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1  Introduction 

The  mixing  of  tracers  on  the  surface  of  a  sphere  has  many  geophysical  applications  from  the 
advection  of  volcanic  ash  to  the  global  redistribution  of  temperature  and  angular  momentum 
by  atmospheric  general  circulation.  Due  to  the  difficulties  of  spherical  geometry,  such  as 
the  poles  and  the  isotropy  between  the  meridional  and  zonal  directions,  previous  work  on 
optimal  mixing  focuses  instead  on  the  d-dimensional  torus  [5,  7,  8].  Here  we  are  motivated 
by  the  geophysical  problem  of  temperature  redistribution  over  the  surface  of  the  sphere 
given  a  heat  source  at  the  equator  and  cooling  at  the  poles.  What  flow  most  efficiently 
mixes  the  temperature  field  and  how  do  the  characteristics  of  these  optimal  flows  compare 
to  similar  flows  on  the  torus? 

Before  describing  the  properties  of  optimal  flows,  we  must  first  define  a  measure  of 
mixing  efficiency  that  will  then  be  used  to  define  a  flow  that  most  optimally  mixes  a  tracer 
field.  One  measure  of  mixing  is  the  variance  of  the  tracer  concentration.  If  we  consider  the 
time  evolution  of  a  tracer  field  with  zero  mean  subject  to  advective  stirring  and  diffusivity, 
the  tracer  concentration  6  obeys  the  advection-diffusion  equation 

r)0 

— -  +  u  •  V0  —  kA 6  =  0  (1) 

at 

where  u  is  the  advective  velocity  (independent  of  the  tracer  concentration  6)  and  where  k, 
is  the  molecular  diffusivity  coefficient.  Multiplying  (1)  by  6  and  averaging  over  the  domain 
yields  the  time  rate  of  change  of  the  variance  of  6  given  by 

^ P -  =  — 2k<|W|2)  (2) 

where  the  integral  over  the  domain  V  is  denoted  by 

(6)  =  yJy0dV.  (3) 

From  (3)  we  see  that  the  variance  of  the  tracer  concentration  ( \6\ 2)  decreases  monotonically 
as  a  function  of  the  molecular  diffusivity  k  and  the  magnitude  of  the  gradients  of  the  tracer 
field.  As  both  k,  and  (|V0|2)  are  positive  definite,  the  variance  of  the  tracer  tends  to  a 
constant  as  the  gradients  of  the  tracer  field  tend  to  zero  in  the  absence  of  sources  or  sinks 
of  the  tracer.  In  other  words,  diffusion  acts  to  destroy  gradients  in  the  tracer  field  and 
homogenize  6. 
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Diffusivity  can  reduce  the  variance  of  a  tracer  field  in  the  absence  of  stirring,  however 
this  process  of  diffusive  mixing  is  quite  inefficient.  We  can  significantly  enhance  the  process 
of  mixing  and  homogenization  by  including  a  stirring  flow  u.  The  mechanical  act  of  stir¬ 
ring  stretches  and  folds  material  lines  of  the  tracer  field,  enhances  the  gradients  of  0,  and 
ultimately  leads  to  a  more  rapid  homogenization  of  the  tracer  field. 

If  there  are  no  source  terms  to  maintain  gradients  in  9  diffusivity  will  lead  to  a  spatially 
homogeneous  tracer  field  over  time.  We  can  obtain  a  nonzero  (|0|2)  as  t  — >•  oo  by  including 
a  spatially  varying  tracer  source  term  in  (1), 


—  +u-  \79-kA9  =  s(x).  (4) 

at 

The  homogenization  of  the  tracer  field  by  mixing  can  then  balance  the  production  of  tracer 
variance  by  the  source  in  steady  state.  It  is  within  this  steady  state  that  we  will  evaluate 
the  effect  of  stirring  on  the  mixing  process. 

As  we  expect  stirring  to  enhance  mixing,  we  anticipate  that  for  the  same  value  of  k 
the  variance  of  a  stirred  tracer  field  will  be  less  than  that  of  an  unstirred  tracer  in  steady 
state.  Thiffeault  et  al.  [7]  measure  the  relative  mixing  enhancement  due  to  stirring  through 
a  measure  given  by  ratio  of  the  unstirred  to  stirred  steady  state  variances, 


•2  (N2) 

m2)‘ 


(5) 


where  9q  satisfies  the  unstirred,  steady  state  diffusion  equation 


— kAOq  s(x) 


(6) 


and  where  6  satisfies  the  stirred,  steady  state  advection-diffusion  equation 


u  •  V0  —  nA9  —  s(x). 


(7) 


A  mixing  efficiency  of  £  =  1  would  imply  that  the  velocity  field  does  not  enhance  mixing 
relative  to  diffusive  processes,  as  would  be  the  case  if  u  were  parallel  to  lines  of  constant 
9  such  that  the  advective  term  u  •  \7 9  of  (4)  vanished.  Mixing  efficiencies  of  £  <  1  will  be 
obtained  if  u  enhances  gradients  of  9  in  steady  state.  Such  small  mixing  efficiencies  are 
counterintuitive  if  we  assume  that  stirring  tends  to  improve  the  homogenization  of  a  tracer 
with  time,  as  is  the  case  for  flows  with  mixing  efficiencies  of  £  >  1. 

By  defining  the  mixing  efficiency  factor  £  we  can  then  create  an  optimization  problem 
to  determine  the  most  efficient  flow  fields  over  the  surface  of  the  sphere  by  maximizing  £. 
Many  approaches  have  been  taken  to  this  problem  in  the  past,  however  all  have  focused 
on  flows  over  the  torus.  There  are  two  primary  methods  of  maximizing  £:  finding  the  flow 
that  most  efficiently  mixes  a  tracer  field  for  a  given  source  function  or,  alternatively,  finding 
the  source  function  that  is  most  efficiently  mixed  by  a  fixed  velocity  field.  Thiffeault  and 
Pavliotis  [8]  consider  the  problem  of  optimizing  the  source  function  for  a  fixed  velocity  field 
and  find  that  these  sources  tend  to  have  maxima  and  minima  located  near  in  regions  of 
strongest  velocity  and  are  aligned  such  that  the  flows  transport  sources  over  sinks. 

Alternatively,  Lin  et  al.  [3]  and  Thiffeault  [6]  maximize  the  mixing  efficiency  for  a  fixed 
source  by  varying  the  velocity  field.  Lin  et  al.  [3]  determine  the  velocity  field  that  most 
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quickly  reduces  a  measure  of  mixing  similar  to  £  in  the  absence  of  diffusivity.  In  contrast 
to  the  transient  problem  considered  by  Lin  et  al.  [3],  Thiffeault  [6]  introduce  a  means  of 
solving  for  the  velocity  fields  that  most  efficiently  mix  a  tracer  concentration  maintained 
by  a  spatially  varying  source  term  and  subject  to  molecular  diffusivity. 

Both  Lin  et  al.  [3]  and  Thiffeault  [6]  find  optimal  velocity  fields  on  the  torus.  Here  we 
will  expand  on  the  approach  of  Thiffeault  [6]  to  consider  optimal  flows  over  the  surface  of 
the  sphere.  We  choose  to  focus  on  three  points  of  comparisons  between  optimal  flows  on 
the  torus  and  the  sphere:  the  issue  of  upper  and  lower  bounds  on  the  mixing  efficiencies, 
the  structure  of  the  most  optimal  flow  configurations,  and  the  scaling  of  the  efficiencies  in 
the  Peclet  number. 

In  Section  2  we  describe  the  numerical  model  we  will  use  to  find  velocity  fields  that 
maximize  £.  We  then  compare  the  analytical  lower  bound  of  the  mixing  efficiencies  on  the 
torus  to  that  of  the  sphere  in  Section  3  and  see  if  the  efficiencies  scale  as  a  function  of  the 
Peclet  number  on  the  sphere  as  they  do  on  the  torus  in  Section  4.  We  will  then  discuss 
the  possibility  of  an  efficient  ‘sweeping’  flow  for  the  sphere  in  Section  5.  In  Section  6  we 
consider  future  directions  of  this  work. 

2  The  Model 

2.1  Euler-Lagrange  Equations 

We  maximize  £  in  a  spherical  geometry  by  varying  the  advective  velocity  field  while  holding 
constant  the  structure  of  the  source  function  and  the  diffusivity.  The  mixing  efficiency  £ 
as  given  in  (5)  has  two  components:  the  variance  of  the  unstirred  steady  state  tracer  field 
(|0o|2)  and  the  variance  of  the  stirred  steady  state  tracer  field  ( \6\ 2).  These  tracer  fields  can 
be  written  as 


6  =  £_1s  0o  =  £o 1s  (8) 

where  C  and  are  the  advection-diffusion  and  diffusion  operators  given  by 

C  =  u  •  V  —  kA  jCq  =  —  ft  A.  (9) 

From  (8)  and  (9),  we  see  that  the  solution  6q  to  the  diffusion-only  differential  equation 
will  be  constant  under  a  change  in  u.  Changing  u  will  only  influence  the  steady  state 
stirred  tracer  field  6  such  that  maximizing  £  is  equivalent  to  minimizing  (|0|2).  In  order  to 
minimize  (|0|2)  with  respect  to  u ,  we  first  define  the  variation  of  (|0|2)  following  Thiffeault 
[6]  ‘  . . 


S(\9\2)=6(\C-18\2)  (10) 

where  we  have  used  6  =  £-1s  from  (8). 

As  we  are  holding  the  source  function  constant,  only  the  integral  operator  C~x  changes 
under  a  change  in  u  such  that 


6(|0|2)  =  2{jC~1s8jC-1s ) 


(11) 
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which  can  be  simplified  using  the  identity  5C  1  =  —C  15CC  1 , 

<5(|0|2)  =  -2  {C-lsC-l8CC~ls).  (12) 

After  integrating  (12)  by  parts,  4(|0|2)  can  be  written  in  terms  of  the  self-adjoint  operator 

A  =  CC\  (13) 

where  H  =  —  u  ■  V  —  kS  is  the  adjoint  of  the  advection-diffusion  operator,  such  that 

<S(|0|2)  =  -2  (A^sSCC-'s).  (14) 

As  the  variation  in  the  advection  diffusion  operator  only  depends  on  the  variation  of  the 
velocity  field,  SC  —  Su  •  V,  we  can  write  the  variation  of  (\9\2)  in  terms  of  Su, 

5(|0|2)  =  -2{A~ls8u  •  V£_1s).  (15) 

The  functional  derivative  of  (|0|2)  with  respect  to  u  then  is 

=  —2  A-'sVjC-'s.  (16) 

ou 

We  could  find  the  extrema  of  (\9\2)  by  finding  the  zeroes  of  (16),  however  this  would 
tell  us  little  about  the  structure  of  the  most  efficient  flows  without  additional  constraints. 
We  place  a  constraint  on  the  energy  and  restrict  our  attention  to  non-divergent  flows  by 
creating  an  extended  functional  of  (|0|2) 

Hu)  =  (I e\ 2>  +  7  «M2>  -  u 2)  +  {vv  ■  u)  (17) 

whose  variation  with  respect  to  u  is 

=  -A-'sVjC-'s  +  ^u-  Vi/  =  0.  (18) 

OU 

Equation  (18)  is  the  Euler-Lagrange  equation  for  the  constrained  optimization  problem 
where  7  and  v  are  the  Lagrange  multipliers  for  the  energy  and  non- divergence  constraints, 
respectively.  Flows  that  satisfy  these  constraints  have  a  kinetic  energy  (I'm]2)  equal  to  a 
constant  value  U 2  and  are  non-divergent  such  that  V  •  u  —  0. 

2.2  Numerical  Model 

There  are  numerous  methods  of  finding  velocity  fields  that  maximize  £  and  minimize  (|#|2). 
Thiffeault  [6],  for  example,  suggests  a  method  to  directly  solve  the  Euler-Lagrange  equations 
given  in  (18)  for  the  velocity  field.  Here  we  consider  an  alternate  means  of  finding  the 
minima  of  ( \9\ 2).  We  perform  a  two-part  optimization  process  based  on  Matlab’s  built-in 
nonlinear  optimization  routine  fmincon.  This  routine,  a  part  of  Matlab’s  Optimization 
Toolbox,  minimizes  (|#|2)  subject  to  the  non-divergence  and  energy  constraints. 

The  fmincon  routine  minimizes  ( \9\ 2)  using  the  iterative  interior  point  algorithm.  At 
each  iteration  of  u  and  C[u],  we  must  solve  the  advection-diffusion  equation  (8)  for  9.  For 
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the  torus  problem,  we  can  solve  (8)  for  6  by  generating  an  invertible,  spectral  advection- 
diffusion  operator  matrix  C[u\.  The  tracer  concentration  6  is  then  obtained  by  solving  the 
system  of  linear  equations  using  this  non-singular  operator.  Due  to  time  limits,  we  were 
unable  to  generate  a  non-singular,  spectral  advection-diffusion  operator  for  the  sphere.  In 
order  to  find  optimal  solutions,  we  instead  solve  (8)  for  6  using  f  solve  where  f  solve  is 
a  built-in  Matlab  routine  that  solves  a  system  of  nonlinear  equations.  Resorting  to  f  solve 
leads  to  an  increased  computational  cost  in  the  spherical  problem  relative  to  the  torus 
problem. 

In  addition  to  solving  for  6  using  f  solve  for  varying  u ,  we  must  consider  that  not  all 
variations  in  u  are  permitted  due  to  constraints.  The  optimization  problem  is  constrained 
in  two  ways:  the  kinetic  energy  of  the  flow  is  held  fixed  at  a  constant  value  and  the  flow 
field  must  be  non-divergent.  The  energy  constraint  limits  the  complexity  and  strength  of 
the  flow  such  that  the  optimal  flow  field  u  in  steady  state  satisfies 

<M2>  =  u2  (19) 


where  U  is  a  constant  equal  to  unity  for  all  optimal  and  sub-optimal  flows  considered  here. 

The  energy  constraint  (19),  along  with  the  coefficient  of  diffusivity  k  and  the  radius  of 
the  sphere  a,  allow  us  to  calculate  the  Peclet  number  for  a  given  flow, 

K 


where  Pe  is  a  nondimensional  parameter  relating  the  relative  magnitude  of  the  advective 
stirring  to  the  molecular  diffusivity  in  the  advection-diffusion  equation  (1)  such  that  mixing 
processes  with  large  Pe  are  advection-dominated.  The  radius  of  the  sphere  is  held  fixed  at 
a  —  1  so  that  the  Peclet  number  of  all  experiments  here  is  a  function  of  the  diffusivity  only. 

The  non-divergence  constraint  V  •  u  can  be  satisfied  by  defining  a  streamfunction  ip 
such  that  the  zonal  and  meridional  velocities  over  the  sphere,  u  and  v  respectively,  satisfy 


1  dp) 
a  dcp 


1  dp) 
a  cos  <p  d\ 


(21) 


where  —i r  <  <p  <  n  is  the  latitude  and  0  <  A  <  2n  is  the  longitude.  It  should  be  mentioned 
that  we  are  using  a  geophysical  coordinate  system  where  <p  —  0  is  at  the  equator  rather  than 
standard  spherical  coordinates  where  the  polar  angle  is  typically  zero  at  the  north  pole. 

As  we  are  interested  in  flows  over  the  surface  of  a  thin  spherical  shell,  it  is  useful  to 
recast  all  operators  and  variables  in  terms  of  spherical  coordinates.  We  use  a  pseudo-spectral 
approach  to  calculating  the  advection-diffusion  operator  C.  Transformations  between  grid 
space  and  spherical  harmonics  in  addition  to  the  calculation  of  derivatives  in  spherical 
coordinates  are  performed  using  Algorithm  888,  a  package  of  Matlab  routines  for  spherical 
harmonic  functions  [1].  Using  Algorithm  888,  we  cast  all  scalars  in  terms  of  their  spherical 
harmonics, 


e{\A)  =  YJYJ  OT( a,m)  (22) 

n=0  m=—n 
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Figure  1:  The  geophysical  source  in  arbitrary  units,  with  heating  (red)  at  the  equator  and 
cooling  (blue)  at  the  poles. 

where  6  is  the  spectral  coefficient  of  0,  n  is  the  total  wavenumber,  m  is  the  zonal  wavenum¬ 
ber,  and  l^Ti(A,/i)  is  the  spherical  harmonic  of  longitude  A  and  of  fi  —  sin</>.  The  spherical 
harmonic  is  given  by 


Y™{  A,M) 


=  e 


imX  Tjm 


Om)- 


where  P™(/i)  are  the  associated  Legendre  functions  defined  by 


(23) 


^m(M) 


1  /  (2n  +  1)  (r 


(24) 


The  wavenumbers  n  and  m  describe  the  structure  of  the  scalar  field  in  latitude  and 
longitude.  In  the  zonal  direction,  m  is  similar  to  the  x-wavenumber  on  the  torus  and  is 
simply  a  superposition  of  Fourier  modes.  The  geometries  of  the  torus  and  the  sphere, 
however,  differ  significantly  in  the  meridional  direction.  Instead  of  a  y— wavenumber  as  on 
the  torus,  spherical  harmonics  are  cast  in  terms  of  a  total  wavenumber  n  where  n  —  m  are 
the  number  of  nodes  in  the  meridional  direction. 

The  geophysical  source  function  we  consider  here  is  given  by  the  n  =  2,  m  —  0  spherical 
harmonic  with  a  spectral  amplitude  of  —1.  The  source  function  can  then  be  written  as 

=  (-l)r2°(A,M)  =  (25) 
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As  n  —  m  =  2  and  m  —  0  for  the  source  function,  s (A,  </>)  has  two  nodes  in  latitude  located 
in  the  midlatitudes  and  has  no  zonal  structure.  The  amplitude  is  chosen  so  that  there  is  a 
source  of  ‘heat’  at  the  equator  and  a  sink  at  the  poles,  as  shown  in  Figure  1. 

2.3  Efficiencies  of  Sub-Optimal  Flows 

The  fmincon  routine  is  an  iterative  process  that  requires  an  initial  velocity  field  sufficiently 
close  to  the  optimal  solution  in  order  to  converge  on  a  minima  in  ( \6\ 2)  with  respect  to  u. 
Fmincon  can  only  identify  local  minima  and,  at  this  time,  we  are  unable  to  state  whether  or 
not  these  optimal  solutions  minimize  ( \6\ 2)  globally.  The  solution  space  is  littered  with  local 
minima  such  that  a  ‘good’  choice  of  an  initial  velocity  field  is  crucial  for  both  convergence 
and  for  the  best  chance  at  finding  a  true  optimal  velocity  field  u. 

In  order  to  determine  what  constitutes  a  ‘good’  initial  guess,  we  begin  our  investiga¬ 
tion  of  optimal  flows  by  determining  the  efficiencies  of  the  sub-optimal  flows  generated  by 
monochromatic  streamfunctions.  Thiffeault  [6]  note  that  the  velocity  fields  that  most  effi¬ 
ciently  mix  a  tracer  in  the  presence  tends  to  be  flows  that  directly  transport  the  tracer  from 
regions  of  excess  6  (‘hot’  regions,  if  we  imagine  the  tracer  is  temperature)  to  sink  (‘cold’) 
regions. 

For  the  geophysical  source  in  Figure  1  and  following  the  observations  of  optimized  source 
functions  in  Thiffeault  and  Pavliotis  [8],  the  most  optimal  flow  would  likely  be  cellular  with 
the  strongest  velocities  over  the  maxima  or  minima  of  the  source  and  nodes  located  near 
the  mean  value  of  s(x).  The  geophysical  source  function  achieves  its  mean  value  in  the 
midlatitudes  such  that  we  would  expect 

1  dip  1  dip 
a  dcpffi  cos  <j>  d\ 

in  these  regions.  Using  (21),  the  streamfunction  ip  of  this  efficient  flow  would  then  likely 
achieve  its  extrema  in  the  midlatitudes  and  have  an  inflection  point  near  the  equator.  In 
terms  of  spherical  harmonics,  a  streamfunction  with  wavenumbers  satisfying  n  —  m  —  1, 
such  as  the  monochromatic  spherical  harmonic  functions  and  Y32,  are  good  choices  for 

Indeed,  ip  oc  Y^  is  the  most  efficient  sub-optimal  streamfunction  at  Pe  =10  over  the 
range  of  monochormatic  streamfunctions  shown  in  Figure  2.  Streamfunctions  proportional 
to  spherical  harmonic  modes  Y™  with  n  —  m  —  1  tend  to  have  the  highest  mixing  efficiencies, 
as  expected,  while  those  streamfunctions  with  n  —  m  —  0  (inter-hemispheric  flows)  are 
slightly  less  efficient.  In  both  the  n  —  m  —  1  and  n  —  m  —  0  flows,  stirring  enhances  mixing 
such  that  £  >  1.  Velocity  fields  described  by  streamfunctions  with  less  zonal  structure,  or 
with  smaller  m,  tend  to  be  more  efficient  than  those  with  larger  values  of  m,  with  an  optima 
in  £  around  m  —  2  for  both  n  —  m  —  1  and  n  —  m  —  0.  Good  initial  conditions  for  ip  in  the 
optimization  algorithm  would  then  be  the  monochromatic  streamfunctions  proportional  to 
y32  and  Yj. 

One  surprising  conclusion  from  Figure  2  is  that  some  flows  with  n  —  m  —  2  are  in  fact 
inefficient.  For  these  flow  fields,  stirring  actually  increases  (|#|2)  relative  to  the  unstirred 
solution  (\6q |2)  resulting  in  £  <  1.  It  is  clear  that  such  flow  fields  would  not  be  good  initial 
guesses  for  the  optimization  problem  as  stirring  acts  to  increase  the  variance  of  the  steady 
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y!)  SPHERE:  Dependence  on  Flow  Structure 
Pe  =  10 


Figure  2:  Sub-optimal  efficiencies  £  for  streamfunctions  of  the  form  ip  oc  stirring 

the  geophysical  source  field  with  n  —  m  —  0  (blue) ,  n  —  m  —  1  (green) ,  and  n  —  m  —  2  (red) 
for  zonal  wavenumbers  m  —  0  to  10  at  Pe  =  10. 

state  tracer  field.  These  inefficient  flows  are  a  unique  feature  of  spherical  geometry  and  of 
our  choice  of  steady,  time-independent  velocity  fields  as  will  be  shown  in  Section  3. 

3  Lower  Bounds  and  Inefficient  Stirring 

In  Section  2.3  we  determined  that  not  all  stirring  by  a  velocity  field  u  enhances  diffusive 
mixing.  We  obtained  a  range  of  sub-optimal  efficiencies  spanning  both  both  positive  and 
negative  values  of  £.  What  is  the  range  of  possible  £  for  the  sphere  and  how  does  this 
compare  to  that  of  the  torus?  We  will  answer  this  in  two  parts,  focusing  on  the  lower 
bounds  and  the  existence  of  inefficient  flows  here  and  on  the  scaling  of  the  most  efficient 
flows  in  Section  4. 

Previous  work  by  Shaw  et  al.  [5]  obtained  analytic  upper  and  lower  bounds  on  £  for  the 
torus  by  optimizing  over  the  tracer  field  6.  For  tracer  fields  stirred  by  a  uniform  velocity 
field,  where  \7u(x,  t)  —  0  at  each  instant  in  time,  and  supplied  by  a  monochromatic  source, 
Shaw  et  al.  [5]  find  that  the  mixing  efficiencies  are  bounded  by 

1  <  £2  <  1  +  Pe2/dL2kg  (27) 

where  d  is  the  dimension  of  the  torus,  L  is  the  domain  size,  and  ks  is  the  wavenumber  of 
the  source.  The  Peclet  number  for  the  torus  is  the  same  as  that  of  the  sphere  given  in  (20) 
with  the  domain  size  L  being  the  characteristic  length  in  the  place  of  the  radius  a. 
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Eff  -2.43M 

BPH£fi|=  (nwnvl««-13 


Max  Streamfunction:  0.69  Contour  Interval:  0.2 
Solution  Found:  Not  Optimized 


(a)  Efficient  Flow,  8  =  2.44 


(b)  Inefficient  Flow,  8  —  0.75 


Figure  3:  Two  of  the  sub-optimal  flows  associated  with  the  parameter  sweep  in  Figure  2 
where  Panel  3(a)  is  associated  with  the  relatively  efficient  flow  of  mode  m  =  2,  n  —  3  and 
where  Panel  3(b)  is  the  inefficient  flow  of  mode  m  =  1,  n  =  3.  The  streamfunction  of  the 
velocity  field  are  plotted  on  top  of  colored  contours  of  the  steady-state  tracer  field  6. 


Shaw  et  al.  [5]  found  that  inefficient  flows  do  not  exist  on  the  torus  for  monochromatic 
sources,  such  as  the  geophysical  source  function  s(x)  oc  (A,  /^)  that  we  consider  here. 
For  such  monochromatic  sources,  stirring  over  the  torus  always  acts  to  decrease  the  ( \9\ 2) 
in  steady  state  relative  to  the  unstirred  solution.  As  we  have  already  seen  from  Figure  2, 
inefficient  flow  not  only  exist  for  the  geophysical  source  on  the  sphere-they  are,  in  fact, 
fairly  common.  What  are  the  properties  of  these  inefficient  flows  that  lead  to  an  increase 
in  the  variance  of  6  relative  to  the  unstirred  solution? 

If  we  compare  the  most  efficient  sub-optimal  flow  of  Figure  2  to  the  most  inefficient  flow, 
we  find  that  the  concentration  of  the  tracer  field  in  steady  state  between  Panel  3(a)  and 
Panel  3(b)  are  quite  different.  The  efficient  flow  has  stagnation  points  centered  near  the 
mean  value  of  the  source  function  (located  in  the  midlatitudes  for  the  geophysical  source) 
and  strong  zonal  velocities  near  the  source  maxima  and  minima.  In  steady  state,  the 
warm  patches  of  tracer  are  advected  off  equator  in  thin  jets  while  cold  regions  are  advected 
equatorwards,  leading  to  efficient  mixing  of  the  tracer  field  relative  to  the  diffusive-only 
solution. 

The  inefficient  flow  in  Panel  3(b)  has  stagnation  points  centered  on  the  equator  at 
the  region  of  maximum  tracer  input  and  strong  flows  parallel  to  lines  of  constant  9.  The 
maximum  ‘heating’  at  the  equator  is  trapped  at  the  stagnation  points  of  the  flow  while  the 
meridional  flow  everywhere  is  relatively  weak.  Without  the  ability  to  remove  excess  tracer 
at  the  equator  and  transport  it  to  the  sink  at  the  poles,  this  flow  tends  to  increase  the 
tracer  variance  in  steady  state  relative  to  the  diffusive- only  solution  and  thus  has  a  mixing 
efficiency  of  £  <  1. 

As  for  inefficient  flows  on  the  torus,  Shaw  et  al.  [5]  identify  that  £  satisfies 

pi  =  (Iflol2)  >  4tt2  (|A~1s|2)  Efc(W27r)~4  l^fcl2  (9R, 

<|0|2>  “  ^  (|V-1s|2)  J2k(Lk/2n)-2\~Sk\2  1  J 
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which  they  then  solve  using  a  variational  problem  in  6  to  maximize  the  stirred,  steady  state 
( \6\ 2)  subject  to  the  steady  state  constraint  ft(|V0|2)  =  (sO). 

The  results  of  Shaw  et  al.  [5]  for  a  monochromatic  source  are  due  in  part  to  the  similarity 
between  the  the  gradient  and  Laplacian  operators  in  Fourier  space,  evident  in  the  similarity 
in  the  numerator  and  denominator  of  the  last  term  in  (28).  For  a  monochromatic  source  of 
wavenumber  kS)  the  Laplacian  and  gradient  operators  are  simply 

V  =  iks  A  =  -k2s  (29) 

so  that  the  Laplacian  operator  is  the  square  of  the  gradient  operator. 

In  terms  of  spherical  harmonics,  however,  the  gradient  operator  is  significantly  different 
from  the  Laplacian.  The  Laplacian  of  the  spherical  harmonic  function  Y™  is  given  by 

A1T(A,m)  =  \,y).  (30) 

While  the  gradient  operator  is  similar  to  that  of  the  torus  in  the  zonal  direction, 

,n)=imY™,  (31) 

the  meridional  derivative  of  Y™( A,  fi)  is  a  recursion  relation 

r\ 

(1  -  i?)  ^ynm(A, y)  =  -nVn+i,mY™+1( A, m)  +  (n  +  iKrrV™  i(A, y)  (32) 

/  2_  2  \  1/2 

where  vn^m  —  (  \n2r^x  )  •  On  f h e  sphere,  the  Laplacian  operator  of  a  monochromatic 

source  function  is  not  equal  to  the  square  of  the  gradient  operator.  The  simplification  of 
the  efficiency  given  in  (28)  from  Shaw  et  al.  [5]  will  not  apply  here.  It  will  thus  be  necessary 
to  find  an  alternative  approach  to  solving  for  the  analytic  lower  bound  on  £  in  spherical 
coordinates. 

4  Upper  Bounds  and  Peclet  Number  Scaling 

The  upper  bound  given  in  (27)  from  Shaw  et  al.  [5]  suggests  that  the  efficiencies  on  the 
torus  scale  linearly  in  the  Peclet  number  for  large  Pe.  Thiffeault  [6]  also  notes  that  for 
optimization  problems  on  the  torus,  the  optimal  efficiencies  scale  quadratically  for  small 
Pe.  Do  these  scalings  in  Peclet  number  hold  for  the  sphere? 

In  order  to  determine  the  efficiencies  of  optimal  flows,  we  must  first  solve  the  optimiza¬ 
tion  problem  given  in  (18).  We  can  use  the  highly  efficient,  sub-optimal  flows  from  Figure  2, 
namely  the  flows  with  i/;  oc  T32(A,  fi)  and  ^  oc  T22(A,  (i)  as  initial  guesses  at  Pe  =  10  to  obtain 
the  two  optimized  flows  shown  in  Figure  4. 

Both  flows  in  Figure  4  have  characteristics  similar  to  the  efficient,  but  sub-optimal  flows 
of  Figure  2.  The  optimal  flows  tend  to  have  stagnation  points  near  the  mean  value  of  the 
source  function  and  both  have  inter-hemispheric  cells  to  transport  the  tracer  from  high 
concentrations  at  the  equator  to  low  concentrations  near  the  poles.  The  efficiency  of  the 
optimized  flow  with  m  —  2  (Panel  4(a))  is  slightly  higher  than  that  of  the  optimized  flow 
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Eff  =  3.2953 

SPHERE:  onemodeNOEN  ex  =  2 


Eff  =  3.1927 

SPHERE:  onemodeNOEN  ex  =  1 


(a)  Optimized  Flow,  £  =  3.30 


Max  Streamfunction:  1 .55  Contour  Interval:  0.3 
Converged:  Optimized 

(b)  Optimized  Flow,  £  =  3.19 
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Figure  4:  Two  optimized  flows  with  initial  guesses  guided  by  the  efficient,  sub-optimal  flows 
y22  and  Y32  shown  in  Figure  2.  Contours  of  the  streamfunction  are  superimposed  on  filled, 
colored  contours  of  the  steady  state  tracer  concentration. 


Figure  5:  Efficiencies  £  —  1  of  optimized  flows  (blue  circles)  and  sub-optimal  flows  (black 
lines)  as  a  function  changing  Peclet  number.  The  efficiencies  of  sub-optimal  flows  are  those 
arising  from  streamfunctions  proportional  to  Y32  (dashed)  and  (solid).  Red  dashed  lines 
are  the  predicted  scaling  for  the  most  efficient  flows  at  low  Pe  (quadratic  scaling)  and  at 
high  Pe  (linear  scaling). 
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Max  Streamf unction:  0.00  Contour  Interval:  0.0 


Figure  6:  The  cellular  source  on  the  sphere  mapped  in  latitude-longitude  coordinates. 

with  m  —  1  (Panel  4(b))  at  £  —  3.30  as  compared  to  £  —  3.19.  Expectedly,  both  flows  have 
efficiencies  greater  than  the  most  efficient  sub-optimal  flow  in  Figure  2. 

While  we  were  able  to  obtain  these  two  optimal  flows  at  Pe  =  10,  it  was  not  without 
difficulty.  As  state  in  Section  2.3,  the  optimization  algorithm  for  the  sphere  is  an  iterative 
process  that  will  converge  on  local  minima,  if  it  converges  at  all.  Unfortunately,  due  to 
time  limitations,  very  few  optimal  solutions  were  found  while  for  varying  Pe.  Those  that 
we  could  identify  are  plotted  in  Figure  5  along  with  the  efficiencies  of  sub-optimal  flows. 
With  the  exception  of  one  point  at  Pe  —  10,  most  all  optimized  solutions  lie  very  near  that 
of  the  initial  guess  (the  mode). 

In  place  of  looking  at  the  optimized  flows,  we  examine  the  scaling  of  the  efficiencies  of 
the  sub-optimal  ‘best  guess’  streamfunctions  with  monochormatic  or  Y%  modal  structure 
instead  of  the  true  optimal  solutions.  From  Figure  5,  we  see  that  the  maximum  efficiencies 
of  the  Y^  and  Y^  modes  appear  to  scale  in  agreement  with  the  observations  of  Thiffeault 
[6]  and  of  Shaw  et  al.  [5]  for  the  torus,  namely  that  the  mixing  efficiencies  scale  as  Pe  for 
small  Pe  and  transition  to  a  Pe 2  scaling  at  high  Pe. 

The  scaling  of  £  with  Pe  even  holds  if  we  move  away  from  the  geophysical  source.  If  we 
consider  a  cellular  source  on  the  sphere  (Figure  6)  and  its  equivalent  on  the  torus,  we  can 
compare  the  scaling  of  the  efficiencies  between  the  two  geometries  directly.  The  efficiencies 
of  sub-optimal  flows  on  the  sphere  (Panel  7(a))  and  both  optimal  and  sub-optimal  flows  on 
the  torus  (Panel  7(b))  are  plotted  against  Pe  in  Figure  7. 

There  are  two  primary  observations  to  note  in  Figure  7:  the  close  proximity  of  optimized 
flows  on  the  torus  to  the  sub-optimal  cellular  flows  and  the  agreement  with  Pe  scaling 
arguments.  For  optimized  flows  on  the  torus,  as  with  the  optimized  flows  on  the  sphere  in 
Figure  5,  the  optimal  solutions  lie  very  close  to  the  ‘best  guess’  sub-optimal  flows.  This 
suggests  that  the  optimized  efficiencies  on  the  sphere  for  the  cellular  source  should  be 
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Figure  7:  Scaling  of  efficiencies  £  —  l  with  respect  to  changing  Pe  for  the  cellular  source 
function  (shown  in  Figure  6  for  the  sphere)  on  the  sphere  (Panel  7(a))  and  for  the  equivalent 
flow  on  the  torus  (Panel  7(b)).  Efficiencies  of  sub-optimal  flows  of  (solid)  and 
(dashed)  over  the  surface  of  the  sphere  are  shown  in  Panel  7(a).  Similarly,  sub-optimal 
flows  on  the  torus  are  shown  for  kx  =  1,  ky  =  1  (solid)  and  kx  =  1,  ky  =  0  (dashed)  in 
Panel  7(b).  Optimized  flows  for  the  torus  are  shown  in  circles.  Quadratic  and  linear  Pe 
slopes  are  red-dashed. 

relatively  close  to  the  sub-optimal  flows  shown  in  Figure  7,  however  this  still  needs  to  be 
shown  in  future  work. 

The  close  agreement  to  the  Pe  scaling  between  the  torus  and  the  sphere  suggests  that, 
unlike  the  lower  bound  on  £,  one  property  of  maximally  efficient  flows  from  the  torus  does 
carry  over  to  the  sphere.  Why,  then,  does  the  efficiency  scale  linearly  in  Pe  at  high  Pe  and 
as  Pe 2  at  small  Pe?  Why  is  this  result  independent  of  geometry? 

At  high  Pe,  we  can  approximate  the  effect  of  stirring  by  and  eddy  diffusivity,  ftejj,  that 
is  significantly  larger  than  the  molecular  diffusivity,  ft.  For  an  advection-dominated  regime, 
an  effective  diffusivity  can  be  generated  by  fteyj  ^  PL,  where  U  is  the  characteristic  length 
of  the  advective  velocity  and  L  is  a  characteristic  mixing  length  [4,  5].  For  the  cellular  flows 
we’ve  considered  here,  where  the  most  efficient  flows  tend  to  be  those  of  the  largest  scale 
permissible,  the  characteristic  length  scale  is  that  of  the  domain  for  the  torus  or  that  of 
radius  for  the  sphere. 

The  stirred  steady  state  ad vect ion-diffusion  equation  (1)  is  then  parameterized  by 

KeffA6  «  s  (33) 

while  the  unstirred,  diffusive  steady  state  problem  is  unchanged 

ftA6>o  =  s  (34) 

so  that  the  solutions  to  the  steady  state  stirred  and  unstirred  systems  are 

e  «  — A-1s  0o  =  -A_1s.  (35) 

Keff  K 
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The  efficiency  £  is  then  approximated  by 


£  = 


QA-^I2) 

(|A-1s|2) 


(36) 


which  scales  as 


£ 


Keff 


UL 


K 


=  Pe 


(37) 


where  we  have  used  Ke$f  ~  UL.  From  (37),  the  mixing  efficiency  £  should  scale  linearly 
with  Pe  so  long  as  Pe  is  sufficiently  large  such  that  the  mixing  is  advection-dominated  and 
the  advection-diffusion  equation  can  be  approximated  by  an  eddy  diffusivity. 

The  diffusion-dominated  small  Pe  regime,  on  the  other  hand,  appears  to  scale  quadrati- 
cally  in  Pe  for  the  cases  we  have  considered  here.  Consider  the  advection-diffusion  operator 
C  given  in  (9)  in  the  limit  of  small  Pe.  By  definition  the  Peclet  number  is  the  scale  of 
the  advection  term  relative  to  the  diffusive  term  in  (1)  and  could  be  considered  a  small 
parameter  in  the  steady-state,  nondimensional  advection-diffusion  equation 


(^Peu  •  V  —  A^j  9  =  s 


(38) 


where  hatted  variables  are  nondimensional.  We  have  assumed  here  that,  to  a  first  approxi¬ 
mation,  the  size  of  the  source  function  scales  as  A 9  in  the  limit  of  small  Pe. 

The  inverse  of  the  nondimensional  operator  £_1  =  ^ Peu  •  V  —  A^  can  then  be  ap¬ 
proximated  by  expanding  about  Pe  <C  1.  Let  A  —  ii  •  V  be  the  advection  operator  and  let 
D  —  —A  be  the  diffusion  operator.  Additionally,  let  e  —  Pe  be  the  small  parameter.  The 
nondimensional  advection-diffusion  operator  (dropping  the  hats)  is  then 


C  =  D  +  eA. 


(39) 


The  inverse  of  (39)  can  then  be  expanded  out  to  O(e) 

C~l  =  D~l  -  eD~lAD~l  +  0(e2) 


such  that  ( \9\ 2)  =  (| C  1s|2)  is 

(\e\2)  =  ±  I  (c~lsfdv 

=  ^J  {D^s-eD-'AD-'s  +  Oie2))2  dV 

Expanding  the  square  of  (41)  yields 

(|0|2)  =  1  J  (D~lsf -2e  (. D~ls )  (yD-lAD-ls)  dV  +  0(e2) 

The  advection  operator  can  be  written  in  terms  of,  =  [h>,  £],  where 

.  ,  dip  dp,  dp  dip  dp  dp 

dx  dy  dx  dy  dy  dx 


(40) 


(41) 


(42) 


(43) 
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on  the  torus,  noting  that  the  non-divergent  advective  velocities  are  defined  using  a  stream- 
function  ip  by  u  —  —  and  v  —  On  the  sphere, 

1  d^\  f  1  \  /  1  d£\  f  1  d'tp 

a  cos  (pdXJ  \ad(p)  \a  cos  <pdXJ  \a  d<p 

for  the  non-divergent  velocities  given  in  (21). 

Replacing  the  advection  operator  in  (42)  with  |/0,£]  yields 

{\9\2)  =  p  I  (D-h)2 -2e(D-1s)(D-1[^,D-1s])dV  +  0(e2)  (45) 

such  that  the  variation  in  ( \6\ 2)  given  a  change  in  the  streamfunction  Sip  is 

<5(|6f)  =  -  y  J  (D-'s)  D~l  [Sip,  D-'s]  dV  +  0(e2)  (46) 

=  -y  J  Sip  [D-'s,  D~2s]  dV  +  0(e2).  (47) 

To  0(e2),  the  variation  of  (\0\2)  is  then 

^p-  =  -2e{[D-1s,D~2s})  +  0(e2).  (48) 

For  any  source  that  satisfies  [D~1s,  D~2s\  =  0,  such  as  monochromatic  sources,  the 
variation  of  (|#|2)  with  respect  to  a  change  in  the  streamfunction  5^>  will  be  zero  at  0(e). 
As  both  the  geophysical  and  the  cellular  sources  considered  here  are  monochromatic,  we  can 
expect  ( \0\ 2)  to  vary  as  e2  —  Pe2,  hence  the  quadratic  scaling  in  the  Peclet  number  at  small 
Pe.  If  we  were  to  instead  look  flows  which  optimized  £  on  the  sphere  with  a  source  function 
satisfying  [P-1s,P-2s]  ^  0,  we  might  expect  linear  scaling  at  small  Pe.  Experiments  so 
far,  however,  have  been  inconclusive.  Finding  sources  where  efficiencies  scale  linearly  in  Pe 
is  left  for  future  work. 

5  Sweeping  Flow 

It  is  interesting  to  note  that  the  flow  appears  to  saturate  in  efficiency  around  Pe  =  102 
in  Figure  5  for  the  geophysical  source.  This  saturation  leads  to  a  crossing  of  the  efficiencies 
of  sub-optimal  flows.  Although  streamfunctions  with  a  modal  structure  of  Y2  are  more 
efficient  mixers  than  those  of  at  small  Pe,  Y^  is  more  efficient  than  the  Y32  mode  at  high 
Pe. 

The  saturation  of  £  at  high  Pe  for  the  Y32  streamfunction  can  be  understood  by  looking 
at  the  steady  state  snapshots  of  0  and  ^  in  Figure  8.  At  Pe  =  103,  the  tracer  field  becomes 
trapped  near  the  stagnation  points  of  u.  For  the  saturated  streamfunction  in  Panel  8(b), 
the  flow  velocity  is  largely  parallel  to  contours  of  constant  0  in  steady  state  such  that  the 
velocity  field  is  no  longer  stirring  the  tracer  field  as,  for  these  flows,  u  •  X7 6  ~  0.  As  the 
advective  term  of  (7)  is  small,  increasing  the  Peclet  number  further  by  either  decreasing  the 
diffusivity  or  increasing  the  strength  of  the  velocity  field  u  will  do  little  to  change  the  mixing 
efficiency,  hence  £  will  not  change  greatly  with  increasing  Pe  as  observed  in  Figure  5. 
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Eti  =  1M.Q591 


SPHCflE-hghPe«eiw-£ 


Max  Streamfunction:  1 .00  Contour  Interval:  0.2 
Solution  Found:  Not  Optimized 


£fl  -  25-B28 
SPHERE  WPESSl  **-  I 


(a)  High  Pe,  S  =  104.06,  Y% 


(b)  High  Pe,  S  =  25.63,  Y£ 


Figure  8:  Flows  given  with  a  streamfunction  structure  of  in  Panel  8(a)  and  in  Panel 
8(b)  at  Pe  —  103  in  the  region  of  efficiency  saturation.  Streamlines  (black)  are  superimposed 
on  colored  contours  of  the  stirred,  steady-state  6. 


We  are  interested  in  finding  the  velocity  field  that  most  efficiently  mixes  the  a  tracer 
over  the  surface  of  the  sphere  for  a  given  source  function,  however  we  have  only  considered 
time  independent  velocity  fields  that  contain  stagnation  points  in  u.  The  presence  of  these 
stagnation  points  limits  the  mixing  efficiency  of  these  steady  flows  at  high  Pe.  We  may  be 
able  to  achieve  higher  £  if  we  allow  the  flow  to  vary  in  time  such  as  to  effectively  remove 
these  stagnation  points  on  average. 

One  example  of  such  an  efficient,  time- varying  flow  for  the  torus  is  found  in  Shaw  et  al. 
[5].  Shaw  et  al.  [5]  are  able  to  saturate  the  upper  bound  of  £  given  in  (27)  using  a  ‘sweeping’ 
flow.  This  flow  is  characterized  by  a  constant  velocity,  (7,  that  alternates  between  a  purely 
zonal  flow  and  a  purely  meridional  flow  over  the  torus  over  long  timescales.  Does  such  a 
sweeping  flow  exist  for  the  sphere? 

A  sweeping  flow  as  described  in  Shaw  et  al.  [5]  for  the  torus  cannot  be  defined  on  the 
sphere  in  the  same  manner  as  it  is  on  the  torus  due  to  the  presence  of  the  poles.  The 
velocity  field  u  must  be  zero  at  the  poles  because  of  artificial  discontinues  in  the  gradient 
operator  in  spherical  coordinates.  If  we  were  to  consider  a  velocity  field  that  was  purely 
zonal  at  one  instant  in  time  it  would  be  require  to  have  the  amplitude  of  u  —  u((/>)  A  to  go 
to  zero  as  (j)  — ►  ±90°.  A  rotation  of  this  flow  by  90°  on  the  sphere  would  not  produce  the 
same  purely  meridional  flow  as  in  the  sweeping  flow  on  the  torus.  Instead,  it  would  produce 
a  flow  with  strong  meridional  velocities  over  the  poles  and  with  two  stagnation  points  at 
the  equator.  With  a  fixed  coordinate  system  this  rotated  flow  violates  the  condition  that 
\u\  —  0  at  the  poles. 

A  sweeping  flow  on  the  sphere  would  then  require  more  than  just  a  rotation  of  the 
velocity  field.  One  method  to  obtain  a  sweeping  flow  on  the  sphere  is  described  by  Finn  [2]. 
Finn  [2]  considers  a  sweeping  flow  produced  by  a  rotation  of  both  the  velocity  field  and  the 
source  function  or,  equivalently,  a  rotation  of  the  coordinate  system. 
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Perhaps  the  simplest  sweeping  flow  over  the  sphere  would  have  an  initial  velocity  field 


u=  f-Ucos^X  (49) 

mixing  a  tracer  supplied  by  the  source  function  given  by 

s(A,  fi)  =  SPx(fi)  (exp  (— iA)  +  exp  (iA))  =  S  (Yj1 (A,  fi)  +  Yf  x( A,  /x))  .  (50) 

This  source  function  is  not  identical  to  the  rotate  geophysical  source.  It  instead  describes  a 
sphere  heated  in  the  western  hemisphere  and  cooled  in  the  eastern  hemisphere.  We  choose 
this  particular  source  instead  of  the  rotated  geophysical  source  in  Finn  [2]  as  simplifies  the 
advective  term  significantly  upon  rotation  of  the  velocity  field. 

We  allow  the  velocity  field  given  in  (49)  to  mix  the  source  from  (50)  to  steady  state 
before  rotating  the  velocity  field  by  90°.  Upon  rotation,  the  velocity  field  will  lie  parallel  to 
lines  of  constant  </>  as  supplied  by  the  source  function  s(A,  /i)  from  (50)  such  that  u  •  V#  =  0. 
As  the  velocity  field  does  not  advect  the  tracer  away  from  the  source,  all  stirring  enhanced 
mixing  occurs  in  the  first  half  of  the  sweeping  flow  when  u  is  given  by  (49). 

The  mixing  efficiency  of  the  sweeping  flow  given  by  (49)  over  the  sphere  can  be  calculated 
by  first  finding  the  steady  state  tracer  concentration  satisfying  the  steady-state  advective 
diffusive  equation  of  (7).  Using  the  gradient  and  Laplacian  formulations  given  in  (31)  and 
(30),  the  tracer  concentration  in  steady  state  satisfies 


E 


n 


(Kv) 


\n=0  m=—n 


+«EE  o 

n=0  m=—n 

=  S{Y1\X^)  +  Y1-1(X,f,)) 


such  that  the  we  can  solve  for  the  9 ™  by  matching  coefficients, 


By  Parsevals’  theorem,  the  variance  of  9  is  then 


(|0|2>=2)r^  El  |C 

n=0  m=— n 


(51) 


(52) 


(53) 


which,  given  the  coefficients  from  (52),  allows  us  to  find  the  steady  state  tracer  variance  of 
the  stirred  problem 


(|6>|2)  =  4t tS2 


3  U2\ 

+  2  If) 


(54) 
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By  taking  U  =  0  in  (54),  we  can  similarly  find  the  steady  state  tracer  variance  of  the 
unstirred  problem 


(N2)  =47 tS2  (4^) 


(55) 


we  can  calculate  the  efficiency  of  the  sweeping  flow  as 

02  (N2)  «4  (A  I 

( \9\ 2)  4 k2  \  a4  2  a2) 


3  o 
1  +  -Pe2 


(56) 


where  Pe  = 

If  we  were  to  compare  the  mixing  efficiency  of  the  sweeping  flow  over  the  sphere  to  the 
upper  bound  for  the  mixing  efficiency  of  a  monochromatic  source  similar  to  that  of  (50) 
on  the  torus  in  (27),  we  would  expect  that  the  most  efficient  sweeping  flow  over  the  torus 
would  have  a  mixing  efficiency  of 

Ztorus  =  1  +  2Pe2  (5^) 

where  the  source  function  on  the  d  —  2  torus  is  described  by 

s(x,  y )  =  St0rus(elx  +  e~lx)  (58) 

in  a  domain  of  length  L  —  1  with  a  total  wavenumber  of  A;2  =  1. 

In  comparing  (56)  to  (57),  we  note  that  the  mixing  efficiency  of  the  sweeping  flow  on 
the  sphere  is  very  close  to  that  of  the  torus.  The  effiency  £ 2  of  the  sweeping  flow  on 
the  sphere  scales  as  | Pe2  as  compared  to  a  scaling  of  \Pe 2  on  the  torus.  For  the  same 
Peclet  number,  the  sweeping  flow  over  the  torus  will  be  slightly  more  efficient  than  the 
corresponding  sweeping  flow  on  the  sphere. 

Additionally,  the  sweeping  flow  over  the  sphere  has  an  efficiency  £  —  yj  1  +  | Pe2  which 
scales  as 


\flPe  for  Pe »  1  (59) 

[  1  +  ^Pe2  for  Pe  «  1 

such  that  £  scales  linearly  in  Pe  at  high  Pe  and  quadratically  in  Pe  at  low  Pe  in  agreement 
with  the  arguments  for  highly  efficient  flows  found  in  Section  4. 

While  Shaw  et  al.  [5]  found  that  the  sweeping  flow  on  the  torus  is  the  most  efficient 
flow  for  a  monochromatic  source  function,  we  are  not  guaranteed  that  the  sweeping  flow 
described  by  (49)  is  the  most  efficient  flow  for  the  the  sphere.  One  particular  issue  with 
any  flow  on  the  sphere  arises  due  to  the  location  of  stagnation  points  in  the  flow. 

As  we  must  have  |tx|  =0  at  the  poles  for  any  rotated  coordinate  system  we  choose,  we 
will  always  have  stagnation  points  in  the  velocity  field  over  the  poles.  It  would  then  always 
be  possible  to  achieve  the  trapping  behavior  described  for  inefficient  flows  in  Section  3  if 
these  polar  stagnation  points  are  located  away  from  the  mean  value  of  the  source  function. 
(One  could  consider  regions  where  the  source  is  equal  to  its  mean  as  already  well-mixed, 
thus  stirring  by  u  could  do  little  to  reduce  the  variance  in  these  regions.)  This  is  one  reason 
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why  we  chose  the  source  given  in  (50)  rather  than  the  rotated  geophysical  source.  The 
rotate  geophysical  source  from  Finn  [2],  given  by 

\/3  1 

s(A,  /x)  =  — — S-PK/k)  (exP  (— 2*A)  +  exp  (2zA))  +  ~^SP2  (m)  ,  (60) 

would  have  a  maxima  over  the  poles.  These  regions  could  never  be  mixed  for  a  fixed  spherical 
coordinate  system,  leading  to  a  saturation  of  £  with  increasing  Pe  similar  to  the  suboptimal 
flows  stirred  by  a  streamfunction  proportional  to  in  Figure  5  and  in  Panel  8(b). 

Finn  [2]  suggests  a  more  efficient  sweeping  flow  could  be  achieved  by  a  careful  rotation  of 
the  coordinate  system  at  each  time  interval.  If  the  poles  of  the  spherical  coordinate  system 
are  always  located  in  regions  where  the  source  function  achieves  its  mean  value  one  could 
bypass  the  trapping  issue  on  the  sphere  and  possibly  find  the  most  efficient,  time-dependent 
flows  on  the  sphere. 


6  Discussion  and  Future  Directions 


In  this  work  we  sought  answers  to  the  two  following  questions:  1)  what  flows  most  efficiently 
mix  a  tracer  field  on  the  surface  of  the  sphere  that  is  supplied  by  a  spatially- varying  source 
function?  and  2)  how  do  the  properties  of  these  flows  compare  to  similar  velocity  fields  on 
the  torus?  We  were  able  to  find  flows  that  maximized  the  mixing  efficiency  given  in  (5) 
using  a  numerical  optimization  scheme.  This  optimization  scheme,  however,  only  optimizes 
the  velocity  field  according  to  one  measure  of  the  mixing  efficiency.  In  determining  the  flows 
that  maximize  other  possible  measures  of  mixing  efficiency,  such  as  the  multi-scale  mixing 
efficiencies  described  in  Shaw  et  al.  [5]  given  by 


<7  2  (IV^ol2) 
p  (|VP0|2) 


(61) 


where  p  =  —1,0,1,  we  could  emphasize  the  homogenization  of  variance  of  6  on  smaller 
(p  —  1)  and  larger  (p  —  —1)  spatial  scales. 

In  changing  geometries  we  find  that  many  of  the  properties  of  optimized  flows  over  the 
torus  do  not  translate  directly  to  that  of  the  sphere.  The  hard  lower  bound  of  1  <  £  of 
Shaw  et  al.  [5],  for  example,  is  violated  for  the  geophysical  source  on  the  sphere  due  to  the 
inescapable  presence  of  stagnation  points  in  u  and  the  trapping  of  the  tracer  concentration 
at  its  source. 

The  issue  of  stagnation  points  in  the  stirring  field  on  the  sphere  also  has  consequences 
for  saturating  an  upper  bound,  should  such  an  upper  bound  be  discovered  in  the  future.  To 
remove  these  stagnation  points  we  must  consider  time-dependent  velocity  fields  similar  to 
the  sweeping  flow  described  by  Shaw  et  al.  [5]  modified  by  the  solid  body  rotation  method 
introduced  by  Finn  [2].  It  may  be  possible  to  determine  the  most  efficient  time-dependent 
flow  for  a  given  source  function  on  the  sphere  by  using  the  method  of  steepest  decent 
described  in  Lin  et  al.  [3]. 

One  similarity  between  efficiencies  of  optimal  (or  sub-optimal,  in  the  case  of  the  sphere) 
flows  for  the  two  geometries  is  that  of  the  scaling  of  £  with  changing  Pe.  For  both  the 
geophysical  and  the  cellular  source  functions,  the  efficiency  £  scales  linearly  with  Pe  at  high 
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Pe  and  quadratically  with  Pe  at  small  Pe.  For  diffusion-dominated  flows  at  small  Peclet 
number  we  expanded  the  advection-diffusion  integral  operator,  £_1  around  the  diffusive 
term  and  determined  that  the  £  ~  Pe 2  scaling  in  this  regime  was  a  consequence  of  our 
choice  of  monochromatic  sources.  It  is  an  open  question  as  to  whether  we  can  obtain  linear 
Peclet  number  scaling  for  non-monochromatic  sources,  such  as  sources  with  off-equatorial 
heating. 

As  we  are  motivated  here  by  the  geophysical  question  of  the  redistribution  of  heat  over 
the  surface  of  the  sphere,  the  optimization  of  £  subject  to  more  complex  source  functions 
are  of  significant  interest.  Some  source  functions  in  particular  are  those  with  off-equatorial 
heating,  as  mentioned,  to  simulate  seasonal  changes  in  solar  insolation.  Additional  source 
functions  we  could  consider  include  those  with  longitudinally  localized  heating  to  simulate 
land-sea  surface  temperature  contrasts. 

Throughout  this  work  we  have  considered  a  velocity  field  that  is  independent  of  tem¬ 
perature,  however  in  a  geophysical  context  these  two  fields  are  not  wholly  independent.  It 
may  be  of  further  interest  to  modify  the  code  to  allow  for  coupling  between  the  velocity 
field  and  the  tracer  field. 

The  energy  constraint  here  follows  from  previous  work  on  the  torus.  However,  this  is 
not  the  only  dynamical  constraint  we  may  wish  to  explore  on  the  sphere.  If  we  were  to 
include  some  form  of  rotation  to  our  model  of  optimized  flows  over  the  sphere,  we  could 
additionally  consider  an  angular  momentum  constraint  in  the  form  of 

M  —  fta2  cos2  (j)  +  ua  cos  (62) 

where  M  is  the  angular  momentum,  ft  is  the  rotation  rate,  a  is  the  radius  of  the  sphere, 
and  u  is  the  zonal  component  of  the  non-divergent  velocity  field  u. 

On  the  Earth,  the  redistribution  of  heat  on  the  sphere  from  large-scale  atmospheric 
stirring  is  not  primarily  achieved  by  the  barotropic  circulations  we  have  considered  here, 
but  instead  by  the  near-zonally  symmetric  circulations  of  the  Hadley  and  Ferrel  cells  in 
height  2  and  latitude  (j)  coordinates.  A  question  of  significant  geophysical  interest  would 
then  be  what  zonally  symmetric  circulation  pattern  in  z-(j)  coordinates  most  efficiently  mixes 
temperature  on  the  sphere?  Furthermore,  can  such  a  flow  be  obtained  using  a  numerical 
optimization  scheme  such  as  the  one  used  here? 
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